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Abstract 


Approximate  expressions  for  geodesic  curves  and  the  geodesic  arc-lengthe 
are  obtained  by  straightforward  methods  which  permit  upper  bounds  of  error  to 
be  established  analytically.  The  errors  are  typically  less  than  1.4  parts  per 
million,  and  even  higher  accuracy  Ib  possible  with  additional  corrections. 

Selected  numerical  examples  are  given,  and  calculated  arc-lengths  are 
compared  with  values  obtained  with  Andoyer's  approximate  formula. 
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Parametric  Formulas  for  Geodesic  Curves 
and  Distances  on  a  Slightly  Oblate  Earth 


1.  INTRODUCTION. 

Hyperbolic  Direction  Finders  at  Very  Low  Frequencies  combine  relatively 
high  accuracy,  *  and  operating  ranges  comparable  to,  and  perhaps  even  larger 
than,  the  earth's  radius.  It  Is  therefore  Important  for  computing  lines  of  posi¬ 
tion  to  Inquire  what  allowance  should  be  made  for  the  fact  that  the  earth  Is  more 
nearly  an  oblate  spheroid  than  a  true  sphere.  If  It  Is  assumed  that  the  "flrst-to- 
arrlve"  components  of  an  electro-magnetic  pulse  travel  from  the  source  to  the 
receiver  by  the  shortest  possible  surface  route,  and  If  It  Is  further  assumed 
that  the  velocity  of  propagation  Is  constant  along  this  path,  then  the  problem 
of  computing  ray  trajectories  and  travel-times  Is  equivalent  to  mathematically 
calculating  geodesics  and  geodesic  arc-lengths. 

As  a  result  of  the  considerable  attention  which  has  been  devoted  to  the  geo¬ 
desic  problem,  several  rather  elegant  approximate  solutions  sre  already  avail- 
2 

able.  In  moBt  cases,  the  accuracy  of  these  solutions  1b  high;  but  the  limits  of 
error  are  somewhat  obscure. 

The  treatment  in  the  following  pages  is  straightforward  to  the  point  of  being 
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elementary;  the  condition  for  minimum  path-length  between  two  arbitrary  pointB 
on  the  spheroid  leads  directly  to  a  problem  in  the  calculus  of  variations,  which 
is  then  converted  to  a  differential  equation  in  spherical-polar  coordinates  0  and  0. 
This  equation  is  easily  solved  to  a  high  degree  of  precision  in  terms  of  ah  "ad¬ 
justed"  co-latitude  angle  ?  ,  thus  giving  the  equation  of  the  geodesics  in  paramet¬ 
ric  form.  Geodesic  arc-lengths  are  also  obtained  in  parametric  form  by  approx¬ 
imate  integration  along  the  geodesic  curve.  The  analytic  approximations  used 
are  shown  to  be  better  than  1.4  parts  in  a  million,  and  are  thus  more  than  ade¬ 
quate  for  purposes  of  very  long  range  radio  location  at  the  present  state-of-the- 
art.  Furthermore,  the  approximations  are  of  sueli  simple  nature  that  if  desired, 
an  even  higher  accuracy  can  be  obtained  in  numerical  cases. 


2.  THE  FIGURE  OF  THE  EARTH,  AND  LATITUDE  CONVERSION  FORMULAS 

Precision  surveying  and  mapping  techniques  refer  all  latitudes  and  longitudes 
to  a  reference  spheroid  which  has  been  chosen  to  approximate  the  figure  of  the 
earth,  but  whose  placement  and  dimensions  are  to  a  certain  extent  arbitrary. 

Such  a  frame  of  reference  constitutes  a  geodetic  "datum".  In  the  United  States 
the  North  American  Datum  of  1927  is  employed  in  modern  work.  In  other  count¬ 
ries  other  datums  are  used,  but  the  problem  of  converting  coordinates  in  one 
datum  to  those  in  another  is  beyond  the  scope  of  these  considerations,  which 
assume  that  the  North  American  Datum  is  extendable  over  the  whole  earth,  and 
represents  the  shape  of  the  earth  with  sufficient  accuracy. 

In  thiB  datum,  the  reference  geoid  is  the  "Clarke  Sphcroid"'of  1866  whose 
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dimensions  are  listed  in  Table  I,  along  with  certain  derived  constants  used  in 
the  analysis  to  follow. 

In  rectangular  coordinates  x,  y,  z  the  equation  of  the  spheroid  is 


»  1 


(1) 


where  the  minor  axis  of  the  spheroid  is  taken  to  coincide  with  the  OZ  axis  of 
coordinates.  In  the  corresponding  polar -spherical  coordinates  (R,  6,  0),  the 
equation  is: 
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sin2e  cos2e  .  i_ 

a5  b2  R2 


(2) 


Hence, 


R 


_ a _ 

V 1+8  cos40 


(3) 


Table  1.  List  of  Constants 


Semi  Major  Asia  (a)  =  6378.  2064  km 

Semi  Minor  Asia  (b)  3  6356.  S838  km 

A 

a  a  -  b  =  21. 6226  km 

A 

a 

.  „  0.00338008 

a 

a  +  b 
2 

»  6367.3951  km 

b 

a 

»  0.99660992 

[b\  2 
\Hl 

«  0.9832  3134 

l!)1 

■  0, 98650849 

a 

b 

»  1.0034016 

( a\2 

\bl 

•  1.  0068148 

8 

s  £  -  1  «  0.0068148 
bz 

U 

a  ,e  2 

s  2  "8  "  0.0033900 

Referring  to  Figure  1,  the  geodetic  latitude  of  the  point  P  is  denoted  <P  . 
(This  symbol  is  not  to  be  confused  with  <p,  the  azimuthal  angle  in  the  polar  co¬ 
ordinate  system. )  The  Y  O  Z  plane  cuts  the  spheroid  in  the  ellipse 


a 


(4) 
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whence  it  follows  that 


dz 

dy 


a2  2 


cotV" 


where  is  the  geocentric  latitude.  But, 


dz 

dy 


tan  ( ff/  2  *  )  * 


cot  <P 


and  hence 


tan  <P  ■  tan  <P 1 


^5  cot  6 


Also, 


cos  9 


t44  cot 2<P  +  b4 
2 


sin  6  * 


Va4  +  b4  tan2  9 


It  follows  from  Eq.  (3)  that 


(5) 


(8) 


(7) 


(8) 


(9) 


(10) 
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Figure  1.  Section  of  Spheroid  by  YZ  Plane 


In  the  analysis  to  follow,  much  use  is  made  of  the  parameter  (J  defined  by 
the  relation: 

sin  6 

sin  ?  s  — - - —  (11) 

Vl  +  8  cos2  6 


From  this, 


cos  £ 


VTT§ 


Vl  +  8 


cpb  & 


cos2  0 


(12) 


and 


tan  6  h 

tan  £  =  ,  =  —  tan  9 

Vi  +8  a 


(13) 


The  inverse  relations  are: 

tan  8  ■  VTTs  tan  5  *  "k  tan  ? 

Vl  +  8 

sin  8  =  —  "tv  sin  ? 

Vl  +  8  sin2  ? 


cos  ? 

COS  8  9  7 - - 

Vl+8  sin2  ? 


and,  using  Eq.  (3), 

R  »  ^=4=  -^1  +  8  sin2  S  =  b-^1  +  8  sin2  ? 


Differentiating  Eq.  (14), 

_ i _ ii  s  yrr? 

cos2 B  d(  cos2  ? 


and  by  Eq.  (16), 

d  8  Vl+8  a/b _ 

d  ?  1  +  8  sin2  ?  1  +  8  sin2  ? 


Finally  it  is  often  desirable  to  convert  directly  from  ?  to  9,  and  vice 
without  using  8  .  In  view  of  Eqs.  (7)  and  (14),  the  needed  relations  are: 


(14) 

(13) 


(16) 


(17) 


(18) 


(19) 


versa. 


tan  9  =  ^  cot  ? 


(20) 
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COS  <f 


sin  </> 


(21) 


(22) 


The  corresponding  inverse  relations  are: 


tan  5  *  $  cot  i 


(23) 


COS  £  * 


1 


(24) 


2  cot  ip 


1 


sin  £ 


V 


h2  2 

1  +  —g  *an 

a 


(23) 


The  two  sets  of  equations  are  symmetric. 

The  qualitative  nature  of  the  angleB  ?  and  ^  -  ip  is  illustrated  in  exaggerated 
scale  in  Figure  2.  It  will  be  noted  that  either 


S  or  (26) 


3.  THE  DIFFERENTIAL  EQUATION  OF  A  GEODESIC,  AND  ITS  SOLUTION 

Consider  a  surface  defined  in  spherical  coordinates  by  the  function 
R  -  R  (  6  ).  Then  an  arbitrary  surface  curve  joining  two  points  and  Pg  on 
the  surface  can  be  defined  by  specifying  0  as  a  function  of  9.  The  curvilinear 
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gK*j 

'm 


I 
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distance  S  between  and  P0  is  then 


S  a 


f2  ne 
p, 


(27) 


where 


1  3  1 


(28) 


It  is  now  a  problem  in  the  Calculus  of  Variations  to  chose  0  ( 6 )  so  that  S  will 
be  a  minimum.  The  desired  function  will  be  one  of  the  extremals  satisfying 
the  Euler  condition,  which,  since  I  does  not  contain  0  explicitly  is  simply: 


d 

cTg 


0 


(29) 
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On  integrating  once, 


In  view  of  Eq.  (28),  this  beoomes: 


where  C  is  a  new  constant.  In  order  to  geometrically  interpret  this  constant, 
a  short  digression  is  now  appropriate. 

In  Figure  3  an  elemental  portion  P^G  of  a  geodesic  from  point  Pj  is  shown. 
The  corresponding  elemental  components  along  the  meridian  and  along  the 
minor  circle  are  designated  P^V  and  P^W  respectively.  The  angle  VPjG  is  desig¬ 
nated  B,  and  is  the  bearing  angle  which  the  geodesic  makes  with  respect  to 
north.  In  the  limit, 

P.  W 

tan  B  -  — —  (33) 

P1V 

Now, 


Pj  W  «  R  sin  6  d  0 


In  this,  the  negative  root  is  chosen  so  that  -  d  B  ,  and  hence  P^V,  are  positive. 
It  follows  that: 
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c _ 

fs in2  ?  -  c2 


(36) 


the  second  step  being  a  consequence  of  Eqs.  (3)  and  (11).  On  re-arranglng  Eq. 
(36),  It  will  be  found  that 


C2  =  sin2  5  sin2  B 


(37) 


In  order  to  avoid  ambiguities  of  sign,  the  convention  will  be  adopted  that,  at 
the  starting  point  of  the  geodesic,  the  positive  roots  of  all  the  radicals  are  to 
be  taken.  It  is  then  clear  from  Eq,  (35)  for  example,  that  the  sign  of  C  is  to  be 
taken  to  be  the  same  as  the  sign  of  tan  B. 

If  is  the  bearing  of  the  geodesic  at  £  =  90°  (that  is,  at  the  equator),  it 
follows  that 

C2  -  sin2  B  (38a) 

At  the  "turn-around"  point  of  the  geodesic  (in  B  '  0) 


,  2  , 
sin  ? 


N,S 


(38b) 


where  g  represents  the  closest  approach  of  the  geodesic  to  the  poles,  and  is 
two  valued.  (For  - 1  <  Beq  <  |  ,  the  angles  B  and  ?N  are  equal.) 

It  1b  now  seen  that  the  constant  C  rather  simply  determines  both  the  bearing 
of  the  geodesic  at  the  equator,  and  the  closest  approach  to  the  poles. 

Returning  now  to  the  differential  Eq.  (32),  the  next  step  is  to  change  from  the 
variable  9  to  the  variable  5  [  see  Eq.  (11)  ]  .  Then 


C  A 


sin  £ 


Vsin2  ? 


d? 


.  (39) 


where 
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£■ 


yr^ 


where 


(45) 


(■  «  !l  [  1  +  S  +  62  a(i  -  a)]'1  [  2  52  (1  -a)  -  63  (2  a2  -  3a+  6) 


-  8 4  a(  y  -  a) 


a(9  -  15a)  S6  a2  ] 


(46) 


in  which  for  brevity, 

a  s  sin2  ?  (47) 


Now,  by  Table  1, 

8  «  0.0068148 


and  hence: 

82*  4.64415  X  10"5 
83  »  3.16490  X  10'7 

84«  2.15881  XIO'9  (48) 

85»  1.4698  X  10”U 
S5»  1.0016X10"13 


Hence  it  is  clear  that  (-  is  very  small,  and  that  the  error  in  substituting  A1  for  A 
will  be  approximately  (■  .  it  may  be  seen  from  Figure  4  that  (-  is  everywhere 
less  than  about  1.4X10  6.  Thus,  to  a  very  high  degree  of  precision  Eq.  (39)  can 
be  written 


1  -  u  sin2  ? 

d$»-C  -  d?  (49) 

sin  ?  V  sin2  ?  -  C2 


Figure  4.  f  as  a  Function  of  | 
o 

It  will  be  noted  from  Eq.  (39)  that  C  ^1.  Inspection  of  the  denominator 

of  Eq.  (40)  shows  that  at  least  for  C  #  0,  the  geodesic  cannot  reach  a  latitude 

2  2 

higher  than  that  corresponding  to  sin  ?  «  C  since  beyond  this  point  the  co- 
efficient  of  d  ?  becomes  imaginary.  (See  also  Eq.  (38b),)  If  $  is  to  continue  to 
increase  "beyond"  the  "turn-around"  point,  the  radical  in  Eq.  (49)  must  change 
sign  at  that  point.  If  d  /  d  ?  and  C  were  plotted  for  a  procession  of  points 
following  along  a  geodesic  which  starts  at  the  point  in  the  Northern  Hemi¬ 
sphere,  and  runs  in  a  more  or  less  North-Easterly  direction,  the  result  would 
be  similar  to  that  depicted  in  Figure  5.  There,  starting  at  P^,  ?  is  decreasing 
and  d  /  d  ?  is  negative  and  decreasing,  going  to  "minus  infinity"  at  the  turn¬ 
around  point  ,  Thereafter  J  increases,  while  d  $  /  d  !J  is  positive  and  de¬ 
creasing  (upper  branch  in  Figure  5).  A  minimum  is  reached  at  5  “  ,  and  then 

d  0  /  d  5  goes  to  "plus  infinity"  at  the  second  turning  point,  T2,  and  so  on. 

If  and  are  the  coordinates  of  P^,  and  0,  {  are  the  coordinates  of  a 
running  point  P  on  the  geodesic,  lying  between  and  T^,  it  follows  from 


and  this  will  be  a  positive  quantity.  If  P  lies  on  the  far  side  of  the  turning  point 
Tj,  the  corresponding  equation  is. 
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where  allowance  was  made  for  the  change  In  sign  of  the  radical.  As  might  be 
expected  it  turns  out  that  the  result  of  integrating  Eq.  (50),  namely 


*  -l 

includes  the  second  case  as  well,  provided  the  functions  sin  are  taken  to  have 
their  "principal  values"  in  the  first  case  (that  is,  ^  <T^),  and  the 

next  larger  values  in  the  region  <  ?  <  Tg,  and  so  on.  The  asterisks  are  used 
as  a  reminder  to  choose  the  correct  range  of  the  quantity  so  marked.  To  avoid 
any  possibility  of  confusion,  reference  may  be  made  to  Table  II. 


$  -I 

Table  II.  Ranges  of  sin 


Range  of  ? 

Range  of  Sin  1 

P1  t0  T1 
?!  Tg 

T2  ‘0  T3 

T3  to  T, 

etc. 

0  to  s/2 

s/2  to  3s/2 

3s/ 2  to  5s/ 2 

5s/2  to  7s/ 2 

etc. 

Notei  Tj,  Tg,  Tg,  etc,  are  the 
successive  "turn-around" 
points  of  the  geodesic. 


Except  in  the  cases  C  »  0,  the  longitude  difference  ^  continually  increases 
as  the  running  point  P  moves  forward  along  the  geodesic,  the  contributions  of  all 
parts  of  the  Integral  are  of  one  sign,  there  being  no  possibility  of  one  part 


cancelling,  or  tending  to  cancel,  another.  Under  this  condition  it  is  justifiable 
to  conclude  that  the  maximum  error  in  0  >  ^  ,  is  less  than  1.4  parts  per 
million.  Thus,  even  if  0  -  0^  is  as  large  as  180°  ,  the  error  in  0  -  0j  Is  less 
than  0,000252°,  which  at  most  would  correspond  to  a  positional  error  along  a 
parallel  of  34  meters.  Of  course,  if  necessary,  the  Integral  in  Eq.  (50),  for 
example,  could  be  broken  up  into  sub-integrals  and  each  corrected  for  the  slight 
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difference  between  A  and  A1-,  thus  reducing  the  error  to  any  desired  limit. 

Eq.  (52)  has  an  important  simplification  when  the  starting  point  P^  of  the 
geodesic  is  on  the  equator,  for  then  the  value  of  the  Indefinite  integral  at  £  »  (Jj 
vanishes,  and 


®ln_1  (^T72  00t?)  "Cu  ®ln"1  (vS) 


(53) 


In  terms  of  the  bearing  angle  Bgq  at  (See  Eq.  (36)  ): 

*  '  *i  8  sin'1  (  tab  Beq  cot  £  )  -  u  sin  Beq  sin'1  ^  jj^  ) 


(54)* 


Eqs.  (52),  (53),  and  (54)  for  the  geodesic  may  now  be  interpreted  geomet¬ 
rically  in  terms  of  a  reference  sphere  (Figure  6)  on  which  the  polar  and  azi¬ 
muthal  coordinates  are  £  and  <j>  respectively.  In  Figure  6,  the  point  A  is  anti¬ 
podal  to  Pj,  and  Fj^A  is  a  great  circle  such  that  the  angle  Z  P^P1  at  P1  is 
B  .  Let  A01  be  the  azimuth  of  P1  and  £  the  polar  angle. 

t 

In  the  spherical  triangle  Z  PjP  ,  the  "sine  law"  gives 


sin  £ 


sin  B  sin  y 

_  e.O  “ 

sin  A  p 


(55) 


y  being  the  angle  Pj  0  P 


1 


The  "coBine  law"  applied  to  the  same  triangle  gives 


sin  y 


cos  £ 
cos  B 


eq 


(56) 


and. 


*.  -1 
■>  sin 

/  cos  £ 

V  cos  B 

\  eq 

(57) 


Eliminating  sin  y  from  Eqs.  (55)  and  (56),  gives  , 


♦Equations  (54)  and  (72)  have  been  programmed  for  an  IBM  1620  computer 
by  the  AFCHL  Technical  Services  Division  under  contract  AF  19(628)-411. 
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sin  A01,  =  tan  Bgq  cot  ? 

and, 

A  ^  !  sin  *  (  tan  B  cot  ?  J 


(68) 


(59) 


The  geodesic,  Eq.  (54),  can  now  be  written 

<p  -  -  u  y  sin  B  (80) 


The  "map"  of  the  geodesic  on  the  reference  sphere  may  therefore  be  constructed 
by  drawing  a  great  circle,  and  then  reducing  the  azimuthal  angle  of  each  point  on 
it  by  uy  sin  Bg  ,  where  y  is  the  angular  measure  of  distance  along  the  great 
circle.  In  Figure  6  the  broken  line  P^PPj  represents  the  geodeBic  derived  from 
the  great  circle. 

Any  geodesic  from  P^  (on  the  equator)  with  0  <  B  <  |  will  intersect  the 
equator  again  at  Pg  which  falls  West  of  the  antipodal  point  A  by  the  small  arc- 
distance 


P„  A  »  IT  u  sin  B  (61) 

ct  cq 

(See  Eq,  (80).)  The  same  geodesic,  if  continued  through  the  southern  hemisphere, 
will  by  symmetry  intersect  the  equator  again  adistance  Sir  u  sin  ,  west  of 

pTa 

P,,  and  so  on.  Unless  the  exact  value  of  -j—  is  a  rational  number,  the  geo¬ 
desic  will  never  close  on  Itself  but  will  continue  to  "creep"  round  and  round  the 
sphere.  Figure  7  is  intended  to  illustrate  qualitatively  the  course  of  geodesic 
from  the  equatorial  point  P^  as  viewed  from  above  the  north  pole  of  an  earth 
with  an  exaggerated  degree  of  flattening.  The  dotted  portions  of  the  curves 
represent  parts  of  the  geodesic  in  the  southern  hemisphere.  The  geodesic  may 
be  visualized  as  the  curve  generated  by  tightly  winding  a  thread  over  the  surface 
of  the  (frictionless)  spheroid,  starting  at  P^  and  passing  through  Pg  . 

In  the  special  case  C  =  0  the  geodesic  starts  out  due  North,  and  it  is  intuitive 
that  it  will  follow  a  meridian.  For  this  case,  Eq.  (53)  correctly  gives 


1 


Figure  6.  The  Reference  Sphere 


Another  special  case  of  interest  is  when  Pj  is  on  the  equator,  and  C — ►  1, 
that  is,  the  initial  bearing  approaches  80°  which  would  take  the  curve  on  an  east¬ 
ward,  equatorial  course.  In  this  case  Eq,  (S3)  is  properly  indeterminate,  since 
on  the  equator  $  may  have  any  value. 

Finally,  there  is  an  interesting  characteristic  of  the  geodesics  in  the  vicin¬ 
ity  of  the  point  A  antipodal  to  the  equatorial  starting  point  Pj.  By  Eq.  (61)  it 

will  be  seen  that  regardless  of  the  equatorial  bearing  angle  B  ,  the  small 

eq 

equatorial  arc  PgA  (See  Figure  7)  can  never  exceed  a  value  of  approximately 
It  u.  This  defines  a  sort  of  limiting  point  such  that  the  angular  distance 
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B'igure  7,  Multi-Turn  Portion  of  a  Geodesic 


L»w  A  w  IT  u 


(62) 


and  there  is  a  similar  point  LE  on  the  east  of  the  antipodal  point.  If  P2  is  pre¬ 
cisely  at  the  antipodal  point,  the  shortest  route  from  P;  to  Pg  is  precisely  over 
a  pole.  If  P2  lies  between  Lw  and  A,  or  A  and  LE,  the  shortest  route  is  neither 
polar  nor  equatorial,  but  is  an  arc  lying  in  either  the  north  or  southern  hemi¬ 
spheres.  However,  if  P2  lies  anywhere  on  the  equator  outside  the  small  arc 

A  ,  the  shortest  route  will  be  precisely  along  the  equator  itself. 

W  E 

4,  ARC  LENGTH  ON  A  GEODESIC 

The  arc  length  along  a  geodesic  between  the  starting  point  Pj_  and  an  arbi¬ 
trary  running  point  P  will  be  found  by  performing  the  integration  indicated  in 
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Eq.  (27),  namely 


where  In  view  of  Eq.  (32),  Eq.  (28)  becomes 


Using  Eq.  (3) 

a!  . _ i _ 

a2  1  +  8  coa2  9 
and  hence  by  Eq.  (11) 

sin2  9  «  sin2  ? 
a2 

On  squaring  Eq.  (17)  and  multiplying  by  Eq.  (IS), 


(63) 


(64) 


(65) 
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On  substituting,  Eq.  (85)  becomes 


Putting  this  in  Eq.  (63)  and  changing  variable  from  6  to  5 


X  A  sin  5 

S  -  -  a  J  ~  "==  d? 

Vein2  C  -  C2 


where 


el 


d  9 

dT 


(86) 


(67) 


is  precisely  the  same  quantity  already  encountered  in  the  previous  section,  (see 
Eq.  (40)  )  where  it  was  shown  that  it  could  be  replaced  by  the  quantity  A'  of  Eq. 
(43)  with  an  error  less  than  1.4  parts  per  million.*  (In  Eq.  (67)  the  minus  sign 

j  2  2 

was  chosen  to  give  a  positive  value  of  S  when  vsin  ?  -  C  is  regarded 
as  starting  out  positive  at  C  =  aB  in  Section  3,  and  again  changing  sign  at  the 
"turn-around"  point.)  Thus 


S 


(  1  -  u  sin2  ?  )  sin  ? 
Vein2  C  -  C2 


dC 


(68) 


On  making  use  of  the  identity: 


*As  mentioned  previously  in  connection  with  the  integral  of  Eq.  (50),  the  range  of 
integration  can  be  broken  into  segments  for  each  of  which  a  correction  can  be 
applied  for  the  slight  difference  between  A  and  A1,  thus  obtaining  an  even  more 
precise  approximation. 


1  -  u  sin2  ?  =  j\  -  |  (  1  +  C2  )J  +  |  £  1  +  C2  -  2  sin2  {"  ,  (69) 


Eq,  (68)  becomes: 


S  a  a  1  -  £  (  1  +  C£  ) 


]  / 


-  Bln  £  d  l 
/sin2  5  -  C2 


.  au  f  L± 

2  J  VsW 

5i 


5-l.:.J.8fn.2?  sin  5  d  ? 


a  {i-|(l  +  C2>}  sV1^ 


-  %  cos  c  Vain2  ?  -  C2 


If  the  starting  point  of  the  geodesic  is  on  the  equator,  ^  ,  and  the 

indefinite  integral  vanishes  at  the  lower  limit,  giving  simply 

S  =  a  |l  -  |(1  +  C2)|  sin'1  J===!  -  fcos?  "Vsin2  ?  -  C2  ( 


In  the  last  two  equations,  the  asterisks  (as  before)  are  reminders  to  change 
from  "principal  values"  to  the  appropriate  branch.  In  this  connection  it  is  noted 
that  in  the  first  Integral  of  Eq.  (70),  the  integrand  is  essentially  positive  since  at 
the  start,  d  ?  is  negative,  and  later  When  the  radical  changes  sign,  so  also  does 
d  £  .  In  the  second  integral,  the  integrand 


1  +  C2  -  2  sin2  ? 


Q  S  12 

1/sin  C  -  C 
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Figure  8.  Qualitative  Behavior  of  Integrand  "Q" 


has  a  behavior  of  the  type  illustrated  qualitatively  in  Figure  8.  Starting  at  Pj 


on  the  equator,  Q  = 


•vrr 


Then  following  along  the  geodesic  Northerly 


and  Easterly,  Q  increases  to  zero  when  sin  IJ  • 


1  +  C 


,  and  then  approaches 


+  «as  the  turning  point  of  the  geodesic  is  reached.  Here  the  radical  changes 
sign,  and  Q  increases  from  -  »,  and  so  on. 


It  follows  from  these  considerations  that 


/. 


Q  d  5  must  have  a  be¬ 


havior  somewhat  as  shown  in  Figure  8.  Since  when  starting  out  along  the  geo¬ 
desic  from  P. ,  d  5  is  negative,  the  integral  at  first  increases,  reaching  a  max- 

1  O 


imum  when  sin  Z  ~ 


1  +  C‘ 


,  and  then  decreases  to  zero  at  the  first  turning 


Figure  9.  Qualitative  Behavior  of  J  ^  Q  d  5 


point  Tj  .  Proceeding  further  along  the  geodesic,  Q  is  negative  and  d  £  is  pos- 

1  +  C2 

itive,  eo  that  the  integral  goes  negative,  reaching  a  minimum  when  sin  ?  -  — ^ - 

It  then  increases  through  zero  (when  the  geodesic  crosses  the  equator)  and  so 

on.  Evidently  then,  the  sign  of  the  quantity  *Vsin2  ?  -  C2  is  to  be  changed  at 
each  reflection  point.  In  Eq.  (72)  the  factor  cos  ?  automatically  takes  care  of 
the  sign  change  at  the  equator. 

In  Eq.  (72),  as  already  known  from  Dec  t' on  3,  the  quantity  sin 

is  simply  the  arc  PjP1  (In  radian  measure),  on  the  reference  Bphere,  Figure  6. 
The  factors  containing  u  are  the  "corrections"  applied  to  this  great  circle 
distance  to  obtain  S  . 
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The  difference  in  the  two  routes  , 

Seq  "  S  *  faU  (1  ‘  Sin  Beq)2'(0^Beq<,T/2  }  (7R> 

can  be  as  large  as  —  a  u  or  about  33.96  km. 

Shortest  Distance  Between  Ordinary  Points  on  the  Equator 

For  ali  pairs  of  points  on  the  equator  excluding  those  considered  abc-e,  it 
is  clear  intuitively  that 

S  =  a  (  A  <j>  )  (79) 


where  (A  $  )  is  the  difference  in  longitude  in  radians.  On  setting 
C  1  1  in  Eq.  (72)  an  indeterminancy  arises,  but  this  can  be  avoided  by  regarding 

the  equatorial  arc  P^2  ^ee  Fig01-6  10)as  the  limiting  case  of  the  geodesic  arc 

P1P21  '  when  th®  point  Pg1  moves  approximately  along  the  meridian  NP,,  to 

approach  P2  .  Thus  for  the  point  Pg1,  the  £  -  value  approaches  ff/2,  and  so 

does  B  ,  so  that  in  the  limit,  Eq.  (54)  for  the  geodesic  curve  becomes: 

eq 


Thus  if  C  and  B  are  varied  so  that 
eq 


k  3  cos  ? 
cos  B 

eq 


1  2 
is  a  constant,  the  point  P2  will  approach  P9.  Then  since  C 


(81) 


sin2  B 


eq 


by 


preceding  page  blank 
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Figure  10.  Limiting  Process  for  Two  Points  on  the  Equator 

Eq.  (38a),  it  follows  from  Eqs,  (72)  and  (80)  that 

S  — »  a  (  1  -  u  )  sin  *  k  =  a  (  A  </>  )  (82) 

and  this  is  the  same  as  Eq.  (79),  as  expected. 


5.  NUMERICAL  EXAMPLES 
5.1.  Meridian  Quadrant 

As  a  first  example,  the  length  of  the  meridian  quadrant  will  be  calculated 

4 

with  the  formula  of  Section  4,  using  the  Bessel  Spheroid,  for  which  the  meridian 
quadrant  is  10,000.8557858  km.  The  semi-major  and  semi-minor  axes  are 


.ufiWifBOi 
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a  »  6377.39715500  km 

b  =  6356.07896325  km 


whence  It  follows  that 

S  »  O.OOB7L9218 


and 


U  «  0.0033428786. 


By  Eq.  (74) 


mq 


la 

~  -j- 


(l-f) 


»  10,000.84923 


(03) 


This  differs  from  the  correct  value  by  only  about  6.5  meters  or  about  6.5  parts 
In  ten-million. 

5.2.  Geodesic  at  B  =  45"  from  Point  on  the  Equator  of  Clarke  Spheroid 


For  this 


C 


Sltl 

eq 


1 

V? 


and  Eq«  (51 )  for  the  geodesic  becomes: 
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-  (*i1  »  sin  *  (cot  £  }  “  sin  ^  ( V"2  cos  £  ) 
*■  V2 


(84) 


5.2.1.  At  the  point  where  the  geodesic  reaches  latitude  ?  =  60'  beyond  the 
first  turning  point,  this  equation  gives 


«  144.7356°  135 

v7 


(85) 


«  144.4120° 


Latitude  £  =  60°  corresponds  to  geodetic  latitude^  =  30.0843°.  By  Eq.  (72) 


S  *  a 


[  < 1  -  •«  u  $  «in'1  (yf) + 1  vn 


0,9974575  X  135  .  a  ftAA,nafr; 

a  L  57:2957795  +  0-°0042373 


2.3508276  a 


»  14,992  •  788  km 


5,2.2.  At  the  point  where  the  geodesic  reaches  latitude  C  =  60°  before  the 
first  turning  point,  Sq.  (84)  gives 


6  -  <b,  *  35.2644°  -  45° 

1  n/2 


=  35.1565° 


(86) 


By  Eq.  (72) 


r'-f; i  . 
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S  *  a  (  1  -  .75  u  )  aln"1  -  ■$  VT^T 

•  a  [^IMr  -  °-00042375] 

•  a  [q.78340129  -  0.00042375] 

«  0.7829775  a 

°  4993.992,  km 


6.2.3.  Length  of  geodesic  arc  on  the  previous  curve,  between  points  where 
?  =  60°  (<p  =  30.0834°  )  is  simply 


14,932.788  -  4,993.BS2  *  9,998.796  km. 

6.  NOTE  ON  ANDOYER'S  FORMULA 

For  calculating  geodesic  arc-lengths,  use  is  often  made  of  a  formula  due  to 
Andoyer,5  which,  in  the  present  notation,  1b 


S  *  a  Cf  +  8  S 


(87) 


where  cr  is  the  great-circle  distance  between  two  points  computed  as  if  the 
earth  were  a  sphere;  that  is. 


COS  <7 

and 


sin  ip ^  sin  ipj  +  cos 


COS  if  g  COB  (  -  03  ) 


(88) 
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8  S 


a  -  b  3  sin  a  -  a 


a  -  b  3  sin  a  +  <s 


(90) 


(91) 


For  two  points  on  the  equator,  ^  ■  <^2  »  0,  and  8S  »  0.  (The  case 

when  o  »  it  ,  for  which  p  is  infinite  will  be  excluded  from  the  present  consider¬ 
ation.)  Then  by  Eq.  (87), 


S  »  au  (92) 

This  is  evidently  the  correct  answer  as  far  as  the  equatorial  route  between 
and  P2  is  concerned.  However,  as  discussed  in  Section  4,  if  the  longitude 
$2  of  P2  is  such  that 

^  +  irj  >  ^[^1  +  v  (1_u)]  (93) 

the  sub-polar  route  is  shorter  than  the  equatorial  route  by  the  amount  [  see 
Eq.  (78)  ]  , 


S  -  S  *  E-iP  (  1  -  sin  B  )2  (94) 

eq  2 

which  can  be  almost  as  much  as  33.96  km.  Thus  m  such  cases,  the  Andoyer 
Formula,  while  still  giving  a  "geodesic  distance",  gives  the  longer  rather  than 
the  shorter  extremal.  In  view  of  these  particular  arcs  of  nearly  180°  for  which 
Andoyer1  s  formula  fails,  it  is  of  interest  to  compare  the  results  of  Andoyer's 
formula  with  those  obtained  by  the  use  of  Eqs.  (71)  or  (72)  in  Examples  5,2, 
5.2.1,  and  5.2.2.  It  will  be  seen  that  the  agreement  is  excellent  in  the  following 
examples. 
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8.1.  Distance  between  (  ij>i  «  Oj  ipj  «  0  >  and  (02  “  144.4120*,  <p2  «  30.0843). 


This  corresponds  to  Example  5.2,  where  it  was  found  that  S  »  14,922.789 
km.  Now, 


oos  a  =  cof)  30.0843°  cos  144.4120°  =  -0.70387249 


a  =  134.7224°  "  2.351349  radians 


3  sin  a  -  3  X  0.7105245  =  2.131573. 


sin2  |  «  -y  (  1  -  cos  <7  )  «  0.85183625 


cos2 1  »  y  (1  +  cosa)  -  0.14818376 


a-b  2.131573  -  2.351349  a-b . 

p-  -gr-  - 0.1T8lg37l -  *  -  V  l.«333189 


a-b  2.131573  +  2.351349  _  a-b  „„„„„„„„„ 
T  - 0.85183025 -  T  5-26285700 


2 

>  S  ■  (p-q)  sin 


Sin2  30.0843°  (1.48333169  +  5.28285700) 


-  21,8226 
8 


(  0. 5012736)*  6.74598889 


-  2.702  825X0.2512752X6.74598869 


«  -  4.5815577 


and 


S  S  a<7  +  8S  *  14,997.39  -  4.58  *  14,992.81  km. 
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This  is  about  0.02  km  longer  than  the  value  of  Example  5.2. 

6.2  Distance  between  ( ^  =  0,  <p t  ■  0)  and  ($2  “  3S*1565°i  f  2  °  30.0843°) 

This  corresponds  to  Example  5.2.1. 

Now, 

cos  a  «  cos  30.0843°  cos  35.1565°  «  0.8652888  X  0.8175822  =>  0.70744472 

a  -  44.9724°  =  0.784918  radians 


12.42630  km 
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and 

S  a  an  t  SS  =  5008.356  -  12,426  •  4993.930  km 

This  answer  is  slightly  less  (0.06  km)  than  that  obtained  in  Example  5.2.1. 

6.3  Distance  between  (0t  «  35.1565°,  <pt  »  30.0843°)  and  ($2»144.412J°, 
<p2  «  30.0843°). 

This  corresponds  to  Example  5.2.2. 

Now, 

cos  C  «  sin2  3Q.0B43°  +  cos2  30.0843°  cos  (144.4120°  -  35.1565°) 

*  0.2512752  +  0.7487247  cos  (  109.2555°  ) 

=  0.2512752  -  0.2469153 

>  0.0043S99 

(7  »  89.7502°  »  1.588436  radians 
3  sin  <7  ■  3  (0.9999905)  «  2.989971g 
sin2 1  ■  j  (l  -cos  (7)  «  0,4978200g 

cos2-|  ■  j  <  1  +  cos  a  )  *  Q,5021799g 

„  h  2.999871.  -  1.566436 
p  ~S~  0.602176-flg 


-  2.854625 
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.  2.999971,  +  1.588436 

a  -  o  5 

8  0.4978200, 


8.172807, 


8  S  *  4  p  sin2  30.0843°  ■  2L|226  x  2.854825  X  0,2512752 


7,754807  km 


S  a  aCT+  8S  =  8891.052  +  7.755  «  9988.807  km 


This  value  Is  only  0.011  km  greater  than  the  value  found  In  Example  5.2.2. 
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